% EW equation model parameters.
% Assume the x domain always ranges from [-pi,pi]. The coefficients are

% Discretization parameters
% Discretization parameters
ModelSettings.SiteNum       = 100;       % Number of sites
ModelSettings.LengthUnit    = 'layer';   % 
ModelSettings.SiteWidth     = 2/sqrt(3); % The round particle

ModelSettings.x_max_EW      = pi;        % The x domain in the EW equation ranges from -pi to pi
ModelSettings.x_min_EW      = -pi;
ModelSettings.dx_EW = (ModelSettings.x_max_EW-ModelSettings.x_min_EW)/(2*ModelSettings.SiteNum);
dx_EW = ModelSettings.dx_EW;
x_EW = ModelSettings.x_min_EW:dx_EW:ModelSettings.x_max_EW;

ModelSettings.x_max         = ModelSettings.SiteNum*ModelSettings.SiteWidth;
ModelSettings.x_min         = 0;
ModelSettings.dx = (ModelSettings.x_max-ModelSettings.x_min)/(2*ModelSettings.SiteNum);
dx = ModelSettings.dx;
x = ModelSettings.x_min:dx:ModelSettings.x_max;

ModelSettings.mode          = Mode;

K_alpha = zeros(ModelSettings.mode,1);
K_beta  = zeros(ModelSettings.mode,1);
for n=1:ModelSettings.mode
    h_sin = sin(n*x_EW);
    dh_sin = zeros(size(h_sin));
    dh_sin(1) = h_sin(1)-h_sin(end);
    dh_sin(2:end) = diff(h_sin);
    K_alpha(n) = sum(dh_sin.^2)/(pi*dx^2);

    h_cos = cos(n*x_EW);
    dh_cos = zeros(size(h_cos));
    dh_cos(1) = h_cos(1)-h_cos(end);
    dh_cos(2:end) = diff(h_cos);
    K_beta(n)  = sum(dh_cos.^2)/(pi*dx^2);    
end
ModelSettings.M2ModeWeighting=[K_alpha,K_beta];

% Parameters about the model
% Assuming the data from KMC is in a domain [-L,L]
ModelSettings.MV      = 'T';

% ModelSettings.FitTo   = 'meanR2';
% ModelSettings.FittingTimeRange = 1000;
% ModelSettings.T       = [300,400,450,500,550,600,620,630,640,650,670,700];
% ModelSettings.W       = 1.0*ones(size(ModelSettings.T));
% ModelSettings.nu      = [2.6068e-003,2.5277e-003,2.4379e-003,1.0669e-003,0.4869e-003,0.2657e-003,9.6627e-005,8.0198e-006,0.8427e-003,6.6350e-003,1.855e-002,4.8547e-002];
% ModelSettings.sigma2  = [0.6307,0.6179,0.6334,0.4737,0.3096,0.1757,0.0291,0.0028,0.0131,0.0395,0.0638,0.0971];
% ModelSettings.Rh      = [1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0];
% ModelSettings.K       = [0.6440 0.6439 0.6520 0.7676 0.9138 0.9823 0.9876 0.9902 0.9909 0.9916 0.9930 0.9939];
% ModelSettings.Tau     = [3.9399 3.9226 4.0073 4.6479 3.5363 2.9746 2.7827 2.6868 2.6295 2.5721 2.4575 2.3667];

% ModelSettings.FitTo   = 'meanR2';
% ModelSettings.FittingTimeRange = 5000;
% ModelSettings.T       = [300,400,500,600,700];
% ModelSettings.W       = 0.2*ones(size(ModelSettings.T));
% ModelSettings.nu      = [3.5407e-004,4.8758e-004,1.9263e-006,6.7600e-004,8.0335e-002];
% ModelSettings.sigma2  = [9.2336e-002,1.2382e-001,6.4740e-003,5.5568e-003,8.9568e-002];
% ModelSettings.Rh      = [0.2 0.2 0.2 0.2 0.2];
% ModelSettings.K       = [0.6440 0.6439 0.7676 0.9823 0.9939];
% ModelSettings.Tau     = [3.9399 3.9226 4.6479 2.9746 2.3667];

% ModelSettings.FitTo   = 'meanM2';
% ModelSettings.FittingTimeRange = 5000;
% ModelSettings.T       = [300,400,500,600,700];
% ModelSettings.W       = 0.2*ones(size(ModelSettings.T));
% ModelSettings.nu      = [5.2895e-006,5.9023e-006,5.4520e-007,5.7229e-005,4.4032e-002];
% ModelSettings.sigma2  = [2.2187e-003,2.4763e-003,2.1068e-004,4.2498e-004,1.3638e-001];
% ModelSettings.Rh      = [0.2 0.2 0.2 0.2 0.2];
% ModelSettings.K       = [0.6440 0.6439 0.7676 0.9823 0.9939];
% ModelSettings.Tau     = [3.9399 3.9226 4.6479 2.9746 2.3667];

ModelSettings.FitTo   = 'meanM2';
ModelSettings.FittingTimeRange = 10000;
ModelSettings.T       = [300,400,500,550,600,620,630,640,650,670,700];
ModelSettings.W       = 1.0*ones(size(ModelSettings.T));
ModelSettings.nu      = [6.5224e-005,8.0291e-005,1.9169e-005,3.2995e-006,1.8432e-007,7.1776e-008,5.6751e-008,4.2134e-004,3.3121e-003,8.1256e-003,1.1854e-002];
ModelSettings.sigma2  = [2.6103e-002,3.1767e-002,9.2163e-003,9.3301e-004,2.7368e-005,4.8094e-006,3.4199e-006,3.2589e-003,2.2196e-002,4.4453e-002,5.1885e-002];
ModelSettings.Rh      = 1.0*ones(size(ModelSettings.T));
ModelSettings.K       = [6.4400e-001  6.4390e-001  7.6760e-001 8.7495e-001 9.8230e-001  9.8462e-001 9.8578e-001  9.8694e-001  9.8810e-001  9.9042e-001 9.9390e-001];
ModelSettings.Tau     = [3.9399 3.9226 4.6479 3.8113 2.9746 2.8530 2.7922 2.7314 2.6707 2.5491 2.3667];

% ModelSettings.FitTo   = 'meanM2+meanR2';
% ModelSettings.FittingTimeRange = 10000;
% ModelSettings.T       = [300,400,500,550,600,620,630,640,650,670,700];
% ModelSettings.W       = 1.0*ones(size(ModelSettings.T));
% ModelSettings.nu      = [3.8983e-003,2.5784e-003,1.4031e-003,3.1008e-004,1.4390e-007,3.3577e-008,3.4466e-008,5.2450e-003,7.5019e-003,2.7020e-002,4.5753e-002];
% ModelSettings.sigma2  = [1.0279e+000,6.6209e-001,5.6588e-001,1.8786e-001,7.4960e-004,9.9617e-005,9.5964e-005,4.6565e-002,4.7125e-002,9.3108e-002,1.0611e-001];
% ModelSettings.Rh      = 1.0*ones(size(ModelSettings.T));
% ModelSettings.K       = [6.4400e-001  6.4390e-001  7.6760e-001 8.7495e-001 9.8230e-001  9.8462e-001 9.8578e-001  9.8694e-001  9.8810e-001  9.9042e-001 9.9390e-001];
% ModelSettings.Tau     = [3.9399 3.9226 4.6479 3.8113 2.9746 2.8530 2.7922 2.7314 2.6707 2.5491 2.3667];

ModelSettings.h = 0;
ModelSettings.rho = 0;
ModelSettings.meanAlpha2 = zeros(ModelSettings.mode,1);
ModelSettings.meanBeta2  = zeros(ModelSettings.mode,1);
ModelSettings.varAlpha2  = zeros(ModelSettings.mode,1);
ModelSettings.varBeta2   = zeros(ModelSettings.mode,1);
ModelSettings.meanR2     = model_CalMeanR2(ModelSettings);
ModelSettings.meanM2     = model_CalMeanM2(ModelSettings);

%%
Settings.D_set     = D_set;
Settings.Fact_D    = Fact_D;
Settings.R2_set    = R2_set;
Settings.Fact_r2   = Fact_R2;
Settings.Ht_set    = Ht_set;
Settings.Fact_H    = Fact_H;
Settings.varR2_set = varR2_set;
Settings.Fact_varR2= Fact_varR2;
Settings.M2_set    = M2_set;
Settings.M2_fact   = M2_fact;
Settings.P         = Step_num;
Settings.m         = Mode;
Settings.dt        = dt;
Settings.lb        = T_LOW;
Settings.ub        = T_HIGH;
Settings.rt_T      = RT_T;
Settings.T_end	   = T_end;
Settings.weight    = [1.0,1.0,1.0,1.0,1.0];
Settings.model     = ModelSettings;
Settings.Input     = {'alpha','beta','h','rho'};
Settings.Output    = {'T(K)'};
Settings.ControllerTypeID = ControllerTypeID;
Settings.FeedThrough = 0;
%%
State.h = 0;
State.rho = 0;
State.meanAlpha2 = zeros(1,Settings.m);
State.meanBeta2 = zeros(1,Settings.m);
State.meanR2 = 0;
State.varAlpha2 = zeros(1,Settings.m);
State.varBeta2  = zeros(1,Settings.m);
State.varR2     = zeros(1,Settings.m);
State.varR2  = 0;
State.meanM2 = 0;
%% Initial guess
T = 500;
